## Mike Ward, and Xun Cao, 29 March  ##
## prediction for gbme model on MID ###
#######################################

out<-read.table("OUT",sep=" ",header=T)
# k scan ll bd1 bd2 bd3 bd4 
#           b0 bs1 bs2 br1 br2
#           s2a sab s2b 
#           s2e rho 
#           s2e1 s2e2 s2e3; s2f1 s2f2 s2f3
n<-dim(out)[1] # of iterations.


burn<-1000 # throw away the first 100,000 iterations
out<-out[-c(1:burn),-c(1,2)]
means<-apply(out,2,mean)

# now get random effects for senders and receivers
d.2000<-dget("data2000")
m<-dim(d.2000[[1]])[1] # of countries in MCMC for this year.

wg.names<-as.character(rownames(d.2000[[3]]))
A<-read.table("A",sep=" ",col.names=wg.names)
B<-read.table("B",sep=" ",col.names=wg.names)
# throw away the first 5,0000 iterations
A<-A[-c(1:burn),]
B<-B[-c(1:burn),]
a<-apply(A,2,mean) # mean random effect for sender
b<-apply(B,2,mean) # mean random effect for receiver
cbind(a,b) # mean sender/receiver effects;
cor.test(a,b)

## now get the mean latent dimensions, z1,z2,z3
V<-scan("V",sep=" ")
dim(V)<-c(3,n,m)
PV<-array(NA,dim=c(m,3,n))
# Put array in correct orientation 
for (j in 1:n)
for (i in 1:m)
    for (k in 1:3){
    PV[i,k,j]<-V[k,j,i]
                }

U<-scan("U",sep=" ")
dim(U)<-c(3,n,m)
PU<-array(NA,dim=c(m,3,n))
# Put array in correct orientation 
for (j in 1:n)
for (i in 1:m)
    for (k in 1:3){
    PU[i,k,j]<-U[k,j,i]
                }

# Find Posterior Mean of v'u
UTV<-matrix(0,m,m)
for (i in 1:(dim(PV)[3]-burn)) { UTV<-UTV+PV[,,i+burn]%*%t(PU[,,i+burn]) }

UTV<-UTV/(dim(PV)[3]-burn)



# Xd has the dyadic variables
Xd<-d.2000[[2]]
bd<-means[2:5]
dyadic.part<-bd[1]*Xd[,,1]+bd[2]*Xd[,,2]+bd[3]*Xd[,,3]+bd[4]*Xd[,,4]

Xr<-d.2000[[3]]
# postscript("pred2000.ps",horizontal = FALSE, family = "Times")
#    postscript("pred2000.ps", 
#                horizontal = FALSE, onefile = FALSE, 
#                family = "Times")
par(mfrow=c(2,1))
par(mar=c(1,3,2,2))
##
c.p<-matrix(NA,nrow=length(wg.names),ncol=length(wg.names))
for (i in 1:length(wg.names)) 
    for (j in 1:length(wg.names))
    {
    c.p[i,j]<-   #  a[i] + b[j] + 
                        means[ 7]*Xr[i,1]+
                        means[ 8]*Xr[i,2]+
                        means[9]*Xr[i,3]+
                        means[10]*Xr[j,1]+
                        means[11]*Xr[j,2]+
                        means[12]*Xr[j,3]
                  #     + UTV[i,j]
    }
theta.hat<-dyadic.part + c.p + means[6]  
quantile(theta.hat, c(.025, .5, .975))
pred<-exp(theta.hat)/(1+exp(theta.hat))
diag(pred)<-0
quantile(pred, c(.025, .5, .975))

Yd<-d.2000[[1]]


plot((as.vector(Yd)), main="Probabilities of MID using only covariates", type="l",ylab="", lwd=6, col="grey", xlab="",xaxt="n",las=1)
lines(as.vector(pred),lwd=6,type="l")



##
c.p<-matrix(NA,nrow=length(wg.names),ncol=length(wg.names))
for (i in 1:length(wg.names)) 
    for (j in 1:length(wg.names))
    {
    c.p[i,j]<-     a[i] + b[j] + 
                        means[ 7]*Xr[i,1]+
                        means[ 8]*Xr[i,2]+
                        means[9]*Xr[i,3]+
                        means[10]*Xr[j,1]+
                        means[11]*Xr[j,2]+
                        means[12]*Xr[j,3]
                        + UTV[i,j]
    }
theta.hat<-dyadic.part + c.p + means[6]  
quantile(theta.hat, c(.025, .5, .975))
pred<-exp(theta.hat)/(1+exp(theta.hat))
diag(pred)<-0
quantile(pred, c(.025, .5, .975))

Yd<-d.2000[[1]]
plot(as.vector(Yd), ylab="",main="Probabilities of MID using complete model", type="l", lwd=6, col="grey", xlab="",xaxt="n",las=1)
lines(as.vector(pred*.98),lwd=6,type="l")
lines(sort(as.vector(pred*.98)),lwd=6,type="l")

#dev.off()





# take a look at how many the model has correctly predicted:
Yd.pred<-pred
thres<-sum(Yd)/(m*(m-1)) # take a threshold: could be any reasonable values between 0 and 1;
Yd.pred[Yd.pred>=thres]<-1; Yd.pred[Yd.pred<thres]<-0; 
add<-Yd+Yd.pred
minus<-Yd-Yd.pred

pred11<-0
c.pred11<-NULL
for (i in 1:dim(Yd)[1]){
for (j in c(1:dim(Yd)[1])[-i]){if (add[i,j]==2){pred11<-pred11+1; c.pred11<-rbind(c.pred11, c(rownames(add)[i],colnames(add)[j]))}
                       }}

pred00<-0
for (i in 1:dim(Yd)[1]){
for (j in c(1:dim(Yd)[1])[-i]){if (add[i,j]==0){pred00<-pred00+1}
                       }}

pred01<-0
c.pred01<-NULL
for (i in 1:dim(Yd)[1]){
for (j in c(1:dim(Yd)[1])[-i]){if (minus[i,j]==-1){pred01<-pred01+1; c.pred01<-rbind(c.pred01, c(rownames(minus)[i],colnames(minus)[j]))}
                       }}
                       
pred10<-0
c.pred10<-NULL
for (i in 1:dim(Yd)[1]){
for (j in c(1:dim(Yd)[1])[-i]){if (minus[i,j]==1){pred10<-pred10+1; c.pred10<-rbind(c.pred10, c(rownames(minus)[i],colnames(minus)[j]))}
                       }}
                       
prediction<-matrix(c(pred11, pred01, pred10, pred00, thres, thres), ncol=3)
colnames(prediction)<-c("predicted 1", "predicted 0", "Chosen Threshold")
rownames(prediction)<-c("actual 1", "actual 0")
print(prediction)                               
xtable(prediction)
sum(Yd); m

c.pred11; # actual MID, and predicted MID;
c.pred01; # no MID, but predicted MID;
c.pred10; # actual MID, and failed to predict.

library(xtable)
prob<-NULL
for (i in 1:dim(c.pred11)[1]){prob<-c(prob, pred[c.pred11[i,1],c.pred11[i,2]])}
Pred11<-cbind(c.pred11, rep(1, dim(c.pred11)[1]), round(prob, digits=5))
Pred11<-as.data.frame(Pred11); colnames(Pred11)<-c("Sender", "Receiver", "Actual MID","PredictedProb")
Pred11<-Pred11[order(Pred11$PredictedProb, decreasing = T),]
xtable(Pred11)

prob<-NULL
for (i in 1:dim(c.pred01)[1]){prob<-c(prob, pred[c.pred01[i,1],c.pred01[i,2]])}
Pred01<-cbind(c.pred01, rep(0, dim(c.pred01)[1]), round(prob, digits=5))
Pred01<-as.data.frame(Pred01); colnames(Pred01)<-c("Sender", "Receiver", "Actual MID","PredictedProb")
Pred01<-Pred01[order(Pred01$PredictedProb, decreasing = T),]
xtable(Pred01)


prob<-NULL
for (i in 1:dim(c.pred10)[1]){prob<-c(prob, pred[c.pred10[i,1],c.pred10[i,2]])}
Pred10<-cbind(c.pred10, rep(1, dim(c.pred10)[1]), round(prob, digits=5))
Pred10<-as.data.frame(Pred10); colnames(Pred10)<-c("Sender", "Receiver", "Actual MID","PredictedProb")
Pred10<-Pred10[order(Pred10$PredictedProb, decreasing = F),]
xtable(Pred10)

Pred<-rbind(Pred11, Pred10)
Pred$PredictedProb<-as.numeric(as.character(Pred$PredictedProb))
Pred<-Pred[order(Pred$PredictedProb, decreasing = F),]

par(mfrow=c(1,1))
plot(rep(c(rep(0,250),1), sum(Yd)), ylab="",main="Actual MID and Predicted Probabiliy, 2000", type="l", lwd=3, col="grey", xlab="",xaxt="n",las=1)
AC<-NULL
for (i in 1:sum(Yd)){AC<-c(AC, c(rep(0,250),Pred$PredictedProb[i]))}
lines(AC*.999,lwd=3,type="l")
